# Import subset from NCCS BMF Data on nonprofits in Syr MSA
load( "../DATA/syr_orgs.Rda" )
# Subset - Extract poverty orgs
pov_orgs <- syr_orgs[which ( syr_orgs$NTEECC == "B61" |
syr_orgs$NTEECC == "B82" |
syr_orgs$NTEECC == "B92" |
syr_orgs$nteeFinal1 == "E" |
syr_orgs$nteeFinal1 == "F" |
syr_orgs$nteeFinal1 == "I" |
syr_orgs$NTEECC == "K30" |
syr_orgs$NTEECC == "K31" |
syr_orgs$NTEECC == "K34" |
syr_orgs$NTEECC == "K35" |
syr_orgs$NTEECC == "K36" |
syr_orgs$nteeFinal1 == "L" |
syr_orgs$nteeFinal1 == "O" |
syr_orgs$nteeFinal1 == "P" |
syr_orgs$nteeFinal1 == "R" |
syr_orgs$nteeFinal1 == "S" |
syr_orgs$nteeFinal1 == "T"), ]print("Don't run me")
# download specific ggmap, register Google key
devtools::install_github("dkahle/ggmap")
register_google(key = 'AIzaSyAMrECdXqpOt443ms6nzh118uUsqC2lE6M',
account_type = "premium", day_limit = 15000)
# subset orgs for gecoding address information
pov_orgs_gps <- subset(pov_orgs, select=c(EIN:INCOME))
pov_orgs_gps$whole_address <- do.call(paste, c(pov_orgs_gps[
c("ADDRESS", "CITY","STATE",
"ZIP5")], sep = ", "))
pov_orgs_gps <- subset(pov_orgs_gps, select=c(whole_address, EIN:INCOME))
for (i in 1:nrow(pov_orgs_gps)) {
latlon = geocode(pov_orgs_gps[i,1])
pov_orgs_gps$lon[i] = as.numeric(latlon[1])
pov_orgs_gps$lat[i] = as.numeric(latlon[2])
}
# Make sure you find a way to then take any PO Boxes and place them in the
# center of centroids, and only the ones that are NAs.head(syr_orgs)load("../DATA/pov_orgs_gps.Rda")
pov_orgs_gps_nona <- na.omit(pov_orgs_gps)
geojson_write( pov_orgs_gps_nona, geometry="point", file="../SHAPEFILES/pov_orgs.geojson" )## <geojson>
## Path: ../SHAPEFILES/pov_orgs.geojson
## From class: data.frame
# I am planning on cutting this section out in favor of doing the buffering a different way in the analysis section
class( pov_orgs_gps_nona )## [1] "data.frame"
coordinates( pov_orgs_gps_nona ) = ~lat+lon
class( pov_orgs_gps_nona )## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
plot( pov_orgs_gps_nona )# Create buffer zone around each point
pov_orgs_gps_nona_buff <- gBuffer( pov_orgs_gps_nona, width=.025, byid=TRUE )
plot( pov_orgs_gps_nona_buff, main="NPO service buffer zone (1.5 miles?)" )class( pov_orgs_gps_nona_buff )## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
Download & clean Census ACS 2011-2015 5-yr data
# Census API
# install.packages("devtools")
# devtools::install_github("hrecht/censusapi")
censuskey <- "f8064a17832b439e690be95ab0e0ef9a78617584"
var.list <- c("NAME",
"B25077_001E", "B25003_001E", "B25003_002E", "B25003_003E",
"B25002_001E", "B25002_002E", "B25002_003E", "B01003_001E",
"B19013_001E", "B17020_002E", "B12006_001E", "B23025_005E", "B15003_017E", "B15003_018E",
"B03001_001E", "B03002_003E", "B03002_004E", "B03002_005E",
"B03002_006E", "B03002_007E", "B03002_008E", "B03002_009E",
"B03001_002E", "B03001_003E",
"GEOID" )
acs_2015 <- getCensus( name="acs5",
vintage = 2015,
key = censuskey,
vars=var.list,
region="tract:*",
regionin="state:36"
)
# subset to Syracuse
# sum(acs_2015$county == "067")
# sum(acs_2015$county == "075")
# sum(acs_2015$county == "053")
acs_2015_syr <- acs_2015[which (acs_2015$county == "067" |
acs_2015$county == "075" |
acs_2015$county == "053"), ]# Race variables are single race, non-hispanic. For example, "white" is really "white non-hispanic." This is true for every race category listed except, obviously, hispanic.
rename( acs_2015_syr,
c("B25077_001E" = "mdn_hous_val",
"B25003_001E" = "tenure_tot",
"B25003_002E" = "tenure_own",
"B25003_003E" = "tenure_rent",
"B25002_001E" = "tot_occup",
"B25002_002E" = "occupied",
"B25002_003E" = "vacant",
"B01003_001E" = "tot_pop",
"B19013_001E" = "mhh_income",
"B17020_002E" = "poverty",
"B12006_001E" = "labor_part",
"B23025_005E" = "unemploy",
"B15003_017E" = "high_sch",
"B15003_018E" = "ged",
"B03001_001E" = "tot_race",
"B03002_003E" = "white",
"B03002_004E" = "black",
"B03002_005E" = "am_ind",
"B03002_006E" = "asian",
"B03002_007E" = "islander",
"B03002_008E" = "other",
"B03002_009E" = "mixed",
"B03001_002E" = "not_hispanic",
"B03001_003E" = "hispanic"
))
acs_2015_syr <- rename.vars(acs_2015_syr,
c("B25077_001E", "B25003_001E", "B25003_002E", "B25003_003E",
"B25002_001E", "B25002_002E", "B25002_003E", "B01003_001E",
"B19013_001E", "B17020_002E", "B12006_001E", "B23025_005E", "B15003_017E", "B15003_018E",
"B03001_001E", "B03002_003E", "B03002_004E", "B03002_005E",
"B03002_006E", "B03002_007E", "B03002_008E", "B03002_009E",
"B03001_002E", "B03001_003E"),
c("mdn_hous_val", "tenure_tot", "tenure_own", "tenure_rent",
"tot_occup", "occupied", "vacant", "tot_pop",
"mhh_income", "poverty", "labor_part", "unemploy",
"high_sch", "ged",
"tot_race", "white", "black", "am_ind",
"asian", "islander", "other", "mixed",
"not_hispanic", "hispanic"))
# Create census race proportion variable
acs_2015_syr$porp_white <- (acs_2015_syr$white) / (acs_2015_syr$tot_race)
names( acs_2015_syr )# Move this chunk to your shapefiles folder to show how you produced the geojson files
# I made another Rmd files with this chunk but also kept the original chunk here. Is that what you would like?
setwd( "../SHAPEFILES/" )
# Census tracts shapefiles of Syracuse MSA, changed to geojson files
Mad_county <- readOGR(dsn="Madison-ct", layer="tl_2010_36053_tract10")## OGR data source with driver: ESRI Shapefile
## Source: "Madison-ct", layer: "tl_2010_36053_tract10"
## with 16 features
## It has 12 fields
## Integer64 fields read as strings: ALAND10 AWATER10
geojson_write( Mad_county, geometry="polygon", file="Mad_county.geojson" )## <geojson>
## Path: Mad_county.geojson
## From class: SpatialPolygonsDataFrame
class(Mad_county)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
plot(Mad_county)
title(main = "Census Tracts in Madison County, NY")names(Mad_county)## [1] "STATEFP10" "COUNTYFP10" "TRACTCE10" "GEOID10" "NAME10"
## [6] "NAMELSAD10" "MTFCC10" "FUNCSTAT10" "ALAND10" "AWATER10"
## [11] "INTPTLAT10" "INTPTLON10"
On_county <- readOGR(dsn="Onondaga-ct", layer="tl_2010_36067_tract10")## OGR data source with driver: ESRI Shapefile
## Source: "Onondaga-ct", layer: "tl_2010_36067_tract10"
## with 140 features
## It has 12 fields
## Integer64 fields read as strings: ALAND10 AWATER10
geojson_write( On_county, geometry="polygon", file="On_county.geojson" )## <geojson>
## Path: On_county.geojson
## From class: SpatialPolygonsDataFrame
class(On_county)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
plot(On_county)
title(main = "Census Tracts in Onondaga County, NY")names(On_county)## [1] "STATEFP10" "COUNTYFP10" "TRACTCE10" "GEOID10" "NAME10"
## [6] "NAMELSAD10" "MTFCC10" "FUNCSTAT10" "ALAND10" "AWATER10"
## [11] "INTPTLAT10" "INTPTLON10"
Osw_county <- readOGR(dsn="Oswego-ct", layer="tl_2010_36075_tract10")## OGR data source with driver: ESRI Shapefile
## Source: "Oswego-ct", layer: "tl_2010_36075_tract10"
## with 30 features
## It has 12 fields
## Integer64 fields read as strings: ALAND10 AWATER10
geojson_write( Osw_county, geometry="polygon", file="Osw_county.geojson" )## <geojson>
## Path: Osw_county.geojson
## From class: SpatialPolygonsDataFrame
class(Osw_county)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
plot(Osw_county)
title(main = "Census Tracts in Oswego County, NY")names(Osw_county)## [1] "STATEFP10" "COUNTYFP10" "TRACTCE10" "GEOID10" "NAME10"
## [6] "NAMELSAD10" "MTFCC10" "FUNCSTAT10" "ALAND10" "AWATER10"
## [11] "INTPTLAT10" "INTPTLON10"
Add Syracuse MSA census tract data to census tract shp files
# drop NAs
Mad_county <- na.omit(Mad_county)
On_county <- na.omit(On_county)
Osw_county <- na.omit(Osw_county)
#Merge spatial files
syr_msa <- union(Mad_county, On_county)
syr_msa <- union(syr_msa, Osw_county)
plot( syr_msa )
# Merge data with MSA shp file
syr_merged <- geo_join(syr_msa, acs_2015_syr, "GEOID10", "GEOID")
plot( syr_merged )
title(main = "Syracuse MSA, NY")# Create geojson file of Syr MSA
geojson_write( syr_merged, geometry="polygon", file="../SHAPEFILES/syr_merged.geojson" )## <geojson>
## Path: ../SHAPEFILES/syr_merged.geojson
## From class: SpatialPolygonsDataFrame
syr_merged_cen = gCentroid(syr_merged,byid=TRUE)
geojson_write( syr_merged_cen, geometry="point", file="../SHAPEFILES/syr_merged_cen.geojson" )## <geojson>
## Path: ../SHAPEFILES/syr_merged_cen.geojson
## From class: SpatialPoints
plot( syr_merged )
points(coordinates(syr_merged),pch=20)
points(syr_merged_cen,pch=20)map <- leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY")
map %>% addProviderTiles(providers$CartoDB)# syr_merged_cen <- geojsonio::geojson_read("syr_merged_cen.geojson", what = "sp")
syr_map <- leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
addPolygons(data = syr_merged,
fillOpacity = 0.3
) %>%
addCircles(lng=pov_orgs_gps_nona$lon, lat=pov_orgs_gps_nona$lat) %>%
addCircles(data = syr_merged_cen)
syr_map %>% addProviderTiles(providers$CartoDB) syr_map <- leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
addPolygons(data = syr_merged,
fillOpacity = 0.3, weight = 2
) %>%
addCircles(lng=pov_orgs_gps_nona$lon, lat=pov_orgs_gps_nona$lat,
color= "red", opacity = 1.0) %>%
addCircles(data = syr_merged_cen, color= "black", opacity = 1.0)
syr_map %>% addProviderTiles(providers$CartoDB) # By race
pal <- colorNumeric(
palette = "Blues",
domain = acs_2015_syr$porp_white )
syr_map <- leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
addPolygons(data = syr_merged, fillOpacity = 0.3, weight = 2,
color = ~pal(acs_2015_syr$porp_white)) %>%
addCircles(lng=pov_orgs_gps_nona$lon, lat=pov_orgs_gps_nona$lat,
color= "red", opacity = 1.0) %>%
addCircles(data = syr_merged_cen, color= "black", opacity = 1.0)
syr_map %>% addProviderTiles(providers$CartoDB)# By Income level
pal2 <- colorNumeric(
palette = "Blues",
domain = acs_2015_syr$mdn_hous_val )
syr_map <- leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
addPolygons(data = syr_merged, fillOpacity = 0.3, weight = 2,
fillColor = ~pal2(acs_2015_syr$mdn_hous_val)) %>%
addCircles(lng=pov_orgs_gps_nona$lon, lat=pov_orgs_gps_nona$lat,
color= "red", opacity = 1.0) %>%
addCircles(data = syr_merged_cen, color= "black", opacity = 1.0)
syr_map %>% addProviderTiles(providers$CartoDB)# URL's below explain how to use functions and arguments in leaflet. It also contains info which notes that the radius of a circle is in meters. Therefore, 1.5 miles converts to 2400 meters.
# https://rstudio.github.io/leaflet/shapes.html
# https://github.com/Robinlovelace/Creating-maps-in-R/blob/master/intro-spatial.Rmd
# Creating a buffer zone, 1.5 miles from every NPO
# I followed the info from this site and tried to apply it to the analysis
# http://walkerke.github.io/2016/12/spatial-pipelines/
# Example using one point:
test_buffer <- c(-76.16674, 43.03231) %>%
st_point() %>%
st_sfc(crs = 4326) %>%
st_transform(32618) %>%
st_buffer(dist = 2400)
test_buffer_map <- test_buffer %>%
st_transform(4326) %>%
as("Spatial")
syr_map %>%
addPolygons(data = test_buffer_map) %>%
addProviderTiles(providers$CartoDB)# Attempt to use every pov_org location, but doesn't work. Coding error.
# pov_orgs_buff <- c(pov_orgs_gps_nona$lon, pov_orgs_gps_nona$lat) %>%
# st_point() %>%
# st_sfc(crs = 4326) %>%
# st_transform(32618) %>%
# st_buffer(dist = 2400)
# pov_orgs_buff_map <- pov_orgs_buff %>%
# st_transform(4326) %>%
# as("Spatial")
# syr_map %>%
# addPolygons(data = pov_orgs_buff_map)
# Alternative way to create buffer for pov_orgs
pov_orgs_buff <- gBuffer(spgeom = pov_orgs_gps_nona, width = 2400)
# proj4string(pov_orgs_buff) <- NA_character_
# proj4string(syr_merged_cen) <- NA_character_
# proj4string(pov_orgs_buff) <- CRS("+proj=longlat +datum=WGS84")
# proj4string(syr_merged_cen) <- CRS("+proj=longlat +datum=WGS84")
# inter <- gIntersection(pov_orgs_gps_nona, syr_merged_cen)
# Chloropleth Map of median housing values
pal2 <- colorNumeric(
palette = "Blues",
domain = acs_2015_syr$mdn_hous_val )
syr_map <- leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addTiles() %>% setView(-76.1474244, 43.0481221, zoom = 9) %>%
addMarkers(lng=-76.1474244, lat=43.0481221, popup="Downtown Syracuse, NY") %>%
addPolygons(data = syr_merged, fillOpacity = 0.3, weight = 2,
fillColor = ~pal2(acs_2015_syr$mdn_hous_val)) %>%
addCircles(lng=pov_orgs_gps_nona$lon, lat=pov_orgs_gps_nona$lat,
color= "red", opacity = .3, radius = 2400 ) %>%
addCircles(data = syr_merged_cen, color= "black", opacity = 1.0)
syr_map %>% addProviderTiles(providers$CartoDB)# Intersection with Census data, capture data
# number of orgs within 1.5 miles of tract's center
# nearest neighbor?
# distance analysis
# Poportion of census tract that are high/low diversity
# number of nonprofits by type in high/low racially diverse areas
# number of nonprofits by type in low/high income census tracts
# Corr analysis between nonprofit access and census tract attributes
# Number of nonprofits location by census tract income or race
# Race
# Income